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Abstract 

Spatially isotropic max-stable processes have been used to model extreme spatial or space-time 
observations. One prominent model is the Brown-Resnick process, which has been successfully 
fitted to time series, spatial data and space-time data. This paper extends the process to possi¬ 
bly anisotropic spatial structures. For regular grid observations we prove strong consistency and 
asymptotic normality of pairwise maximum likelihood estimates for fixed and increasing spatial 
domain, when the number of observations in time tends to infinity. We also present a statistical 
test for isotropy versus anisotropy. We apply our test to precipitation data in Florida, and present 
some diagnostic tools for model assessment. Finally, we present a method to predict conditional 
probability fields and apply it to the data. 
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Keywords: anisotropic space-time process; Brown-Resnick space-time process; hypothesis test for spa¬ 
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1 Introduction 

Max-stable processes, such as the Brown-Resnick process, have been successfully fitted to time series, 
spatial and recently to space-time data. Methods for inference include pairwise likelihood based on the 
bivariate density of the models (cf. Padoan et al. [27]), censored likelihood (cf. Wadsworth and Tawn 
[32] j or threshold-based approaches (cf. Engelke et al. [E]). In Davis et al. |8| a spatially isotropic 
Brown-Resnick space-time process is suggested and applied to precipitation data. Pairwise maximum 
likelihood estimates are shown to be strongly consistent and asymptotically normal, provided the 
domain of observations increases jointly in space and time. Their approach is restricted to isotropic 
spatial dependence. 

In the present paper we generalise the Brown-Resnick model to allow anisotropy in space. The 
new model allows for different extremal behaviour along orthogonal spatial directions. Anisotropy is 
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often observed on Earth, for example in Middle Europe with its westerly winds or near the equator 
where trade winds involve predominant easterlies. All dependence parameters are summarised in the 
semivariogram of an underlying Gaussian space-time process. This semivariogram then defines the 
dependence structure of the max-stable process and, as a consequence, the tail dependence coefficient 
between two process values evaluated at two location and two time points. 

Eurthermore, since in real world applications, observations are often recorded over a large number 
of time points, but only at a comparably small number of spatial locations, we consider both a fixed 
and increasing spatial domain in combination with an increasing temporal domain. Eor both settings, 
fixed and increasing spatial domain, we prove strong consistency and asymptotic normality of the 
pairwise maximum likelihood estimates in the anisotropic model based on regular grid observations. 
This requires in particular to prove space-time and temporal mixing conditions in both settings for 
the anisotropic model. 

We also provide tests for isotropy versus anisotropy again in both settings, which are designed 
for the new model. The asymptotic normality of the parameter estimates determines in principle the 
rejection areas of the test. However, the covariance matrices of the normal limit laws are not available 
in closed form. We formulate a subsampling procedure in the terminology of the Brown-Resnick space- 
time process and prove its convergence for fixed and increasing spatial domain. 

We conclude with an analysis of space-time block maxima of radar rainfall measurements in Elorida. 
Eirstly, we present a simple procedure to test whether they originate from a max-stable process. As this 
cannot be rejected, we fit the Brown-Resnick space-time model to the data, using pairwise maximum 
likelihood estimation. Subsequently we apply the new isotropy test. Both the estimation and the test 
are based on the setting of a fixed spatial domain and increasing time series. In particular, since the 
Brown-Resnick space-time process satishes the strong mixing conditions for increasing spatial and 
time domain as well as for fixed spatial and increasing time domain, the estimation and test procedure 
are independent of the specific setting: it works in both settings in exactly the same way, taking the 
different asymptotic covariance matrices into account. Einally, we assess the goodness of fit of the 
estimated model by a simulation diagnostics based on a large number of i.i.d. simulated anisotropic 
Brown-Resnick space-time processes. As a result, there is no statistical significance that the anisotropic 
Brown-Resnick space-time process with the fitted parameters should be rejected. 

Our paper is organised as follows. In Section we present the Brown-Resnick space-time model, 
which allows for anisotropic effects in space, and various dependence measures, including the param- 
eterised dependence function. In Section we compute the pairwise maximum likelihood estimates 
for the new model and prove their strong consistency and asymptotic normality for both settings, 
fixed and increasing spatial domain. Section presents hypothesis tests for spatial isotropy and de¬ 
rives rejection areas based on a subsampling procedure. A data analysis is performed in Section 
with focus on model assessment. The isotropy test rejects spatial isotropy for these data in favour 
of our new anisotropic model. Based on two other test procedures, we conclude that the anisotropic 
Brown-Resnick space-time process with the given dependence parameters is an appropriate model 
for the block-maxima data. We conclude by predicting conditional probability fields, which give the 
probability of a high value (for example of the amount of precipitation) at some space-time location 
given a high value at some other location. 
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2 Spatially anisotropic Brown-Resnick processes 

Throughout the paper we consider a stationary Brown-Resnick space-time process with representation 

oo 

vis, i) = V } , (s, i) G X [0, oo), (2.1) 

i=i 

where {^j : j G N} are points of a Poisson process on [0, oo) with intensity the dependence 

function 5 is nonnegative and conditionally negative definite and {Wj{s,f) : s G G [0,oo)} are 
independent replicates of a Gaussian process 

{W{s,t) : s G G [0, oo)} with stationary increments, VP(0,0) = 0, E[lP(s,t)] = 0 and covariance 
function 




Representation (2.1) goes back to de Haan m and Gine et al. [I8j Brown-Resnick processes have 
been studied by Brown and Resnick [3] in a time series context, as a spatial model by Kabluchko 
et al. [23], and in a space-time setting by Davis et al. |7j and Huser and Davison [T9|. The univariate 
margins of the process rj follow standard Frechet distributions. 


There are various quantities to describe the dependence in (2.1): 


In geostatistics, the dependence function 6 is termed the semivariogram of the process {W(s, t)}: 
For G X [0 , oo), it holds that 


Yar[Wis^^\t^^^) - Wis^^\t^^^)] = 26{s^^^ - t(i) - t(2)). 


• For G and u G M, the tail dependence coefficient x(ft., u) is given by (cf. Kabluchko et al. 
[23] . Remark 25 or Davis et al. [7], Section 3) 

Xih,u) := l^Fir]{s^^\t^^'l) > y \ r]{s^‘^\t^^'l) > y) = 2^1 - $^ , (2.2) 

where h = — s^'^\ u = and ‘h denotes the standard normal distribution function. 

• For D = {(s(i), ..., (sd^l), fd^D)} and y = (yi,..., y\D\) > 0 the finite-dimensional margins 

are given by 

< yi,y{s^‘^\r^^) < y2, ■ ■ .,ry(s(l^l\< y|£,|) = (2.3) 

Here Vd denotes the exponent measure, which is homogeneous of order -1. 

• The extremal coefficient for any finite set D C x [0, oo) is defined through 

F{r]{s^^\t^^'l) < y, r/(s^^\ t^^^) < y,..., < y) = y > 0; 

i.e., fo = Rd( 1, ■ ■ •, !)■ If \D\ = 2, then (cf. Beirlant et al. [T|, Section 9.5.1) 

X(s(^) - s^‘^\t^^'> - t(^)) = 2-^d- 
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In this paper we assume the dependence function 5 to be given for spatial lag h and time lag u by 


6{h,u) = Y,Cj\hj 

i=i 


^^+Cd+i\u\ 


Oid+1 


{h,u) = {hi,...,hd,u) G 


(2.4) 


with parameters Cj > 0 and aj G (0, 2] for j = 1,..., d + 1. 


Model (2.4) allows for different rates of decay of extreme dependence in different directions. This 


particularly holds along the axes of a d-dimensional spatial grid, but also for other directions. For 
example in the case d = 2, the decreases of dependence along the directions (1,2) and (2,1) differ. 


Model (2.4) can be generalised by a simple rotation to a setting, where not necessarily the axes, but 


other principal orthogonal directions play the major role. The rotation angle then needs to be estimated 
together with the other model parameters. A similar approach has been applied to introduce geometric 
or zonal anisotropy into a spatial isotropic model (see e.g. Blanche! and Davison [2], Section 4.2, or 


Engelke et al. m, Section 5.2). For a justification of model \2A\ see Buhl [H], Sections 4.1 and 4.2. 
There it is shown that Brown-Resnick processes with this dependence function arise as limits of 
appropriately rescaled maxima of Gaussian processes with a large variety of correlation functions. 


3 Pairwise maximum likelihood estimation 


We extend the pairwise maximum likelihood procedure described in Davis et al. [8] for spatially 
isotropic space-time Brown-Resnick processes to the anisotropic case. We focus on the difference 
introduced by the spatial anisotropy and refer to the corresponding formulas in Davis et al. [8], where 
also a short introduction to composite likelihood estimation and further references can be found. 

The pairwise likelihood function uses the bivariate distribution function of 
+ h,t + u)) = {T]{ 0 , 0 ),r]{h,u)) (equal in distribution by stationarity) for h G and 

u gM., which is given as 


G{yi,y2) = exp{-l/(yi,y2)}, yi,y2 > 0, 


(3.1) 


where the exponent measure V 


Vd for D = has the representation 


y{yi,y2) 

^ log(y2/yi) , , IS{h,u)\ , 1 log{yi/y2) , , l6{h,u)\ 


(3.2) 


which is a particular form of Eq. (2.7) in Hiisler and Reiss 
by (2.4). For a derivation of (3.2) see for instance Oesting 


. The dependence function <5 is given 
, Satz und Definition 2.4. 

From this we can calculate the pairwise density g{yi,y 2 ) = ge{yi-,y 2 ) of G by differentiation. The 
parameter vector 6 = (Ci,..., Cd+i-, ai,..., ctd+i) lies in the parameter space 


0 := {(Cl,..., Crf+i, oi,.. 


.,ad+i) : Gj G (0,oo),Q;j G (0,2],j = . 


We focus on data on a regular spatial grid and at equidistant time points. More precisely, we 
assume that the spatial observations he on a regular d-dimensional lattice, 

Sm = {^ = ('Sij • • ■, Sd) G {1,..., M}*^} 
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for M G N, and that the time points are given by the set Tt = {I, ■ ■ ■ ,T} for T G N. 

For the computation of the pairwise likelihood it is common not to include observations on all 
available space-time pairs, but only on those that lie within some prespecified spatio-temporal distance. 
This is motivated by the fact that pairs which lie sufficiently far apart in a space-time sense have 
little influence on the dependence parameters, see Nott and Ryden [25], Section 2.1. To express this 
notationally, we take inspiration from that paper and use a design mask adapted to the anisotropic 
setting, 

Tir ■■= {h = {hi,...,hd) €NQ:h<r}, r = (n,..., r^) G Nq. (3.3) 

We are now ready to define the pairwise log-likelihood function and the resulting estimate. 

Definition 1 (Pairwise likelihood estimate). The pairwise log-likelihood function based on space-time 
pairs, whose maximum spatial lag is r G Nq and maximum time lag is p € No, such that {r,p) / (0, 0), 
is defined as 

T p 

PL(^’^)(0) : = EE E E l{(h,«) 7 ^(o,o)} log {ge {rjis, t),r]{s + h,t + u))} 

1 hGTir u—0 
s+hGSj\^j t-\-u<T 

T 

= J2^eis,t;r,p )0 G 0, (3.4) 

where 

p 

qe{s,t;r,p) := EE l{(Ai,«)^(o,o)} log {ge {v{s, t), ri{s + h,t-t- u))} (3.5) 

h^Tir U=0 

and 

T p 

7^(M,r)(0) _ 

E EE Ell s+h^Sm or t+u>T} log {ge {v{s, t),ri{s + h,t + u))} 

h^^htr 'R—0 
p 

= EE E log {ge {v{s, t),r]{s + h,t + u))} , (3.6) 

hG'Hr- ri=0 (s,bs6M,T(^,ii) 

with 


GM,T{h, u) ■= {(s, t) G Sm xTt '■ s -\-h ^ Sm ot t + u>T} . 
for {h,u) G The pairwise maximum likelihood estimate (PMLE) is given by 


6 = argmaxPLl^’^l(0). 
0ee 


(3.7) 


(3.8) 


We derive the asymptotic properties of the PMLE for two scenarios. The first one is based on 
regularly spaced observations with an increasing spatio-temporal domain. For this scenario we follow 
the proofs in Davis et al. |8] and show that the properties of strong consistency and asymptotic 


normality also hold if the dependence structure 5 allows for spatially anisotropic effects as in (2.4). In 


the second scenario, the observations are taken from a fixed spatial domain and an increasing temporal 
domain. 
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3.1 Increasing spatio-temporal domain 
Lemma 1. For (h, u) G Hr x {0, ... ,p}, it holds that 

\gM,T{h,u)\ < K2{M'^-^T + 

where K 2 is a constant independent of M and T. 

Proof. The number of space-time points within the space-time observation area, from which some grid 

d 

point outside the observation area is within a lag {h, u) G Hr x {1,..., p}, is bounded by M^~^T ^ rj + 

i=i 

M^^p. Thus we obtain 

d 

\gM,T{h, w)l < M^-^T rj + M'^p < K2{M'^-^T + M'^), 

i=i 

where K 2 := max{ ^ constant independent of M and T. □ 

Theorem 1 (Strong consistency for large M and T). Let {p{s,t) ■ s G G [0,oo)} be a Brown- 


Resnick process as in (2.1) with dependence structure 

d 


s{h,u) = YCj\hjr^+Cd+i\ur+^, {h,u) 


G 


pd-\-l 




where 0 < aj < 2 and Cj > 0 for j = 1,... ,d + 1. Denote the parameter vector by 

0 {C\^ • ■ • ? 0^1? • • ■ ; 

Assume that the true parameter vector 0* lies in a compact set 

0* C {(Cl,..., Cd+i, ai,..., Orf+i) : Cj G (0, 00 ), Oj G (0, 2], j = 1,..., d -|- 1} . 
Suppose that the following identifiability condition holds for all {s,t) G Sm x Tt-' 

e = e ^ 

g 0 {p{s,t),'n{^ + I^H + u)) = 99{'n{^H),v{^ + f^H + u)), h^Hr, 0<u<p. 
Then, the PMLE 


(3.9) 


(3.10) 


e 


(M,T) 


= aigmax PL^^C) f0\ 

eee* 


'^{M,T) a.s. 

6 0 as M,T ^ 00 . 


is strongly consistent: 


Proof. The proof uses the method of Wald |33] . One aim is to show that for some chosen maximum 
space-time lag {r,p) G \ {0} and 6 G 0*, 


M<^T 


L^PL(".^)(e) 

M^( ^■PL(^) :=E[g0(l,l;r,p)] 




as M, T —7- 00 . This is done by verifying the following two limit results: Uniformly on 0*, 
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(A) E ^g 0 (s,i;r,p) PL(0) as M,r^ oo, 

sGSm t=l 

(B) "4- 0 as M,T ^ oo. 

Furthermore, we need to show: 


(C) The limit function PL(0) is uniquely maximized at the true parameter vector 6* G 0*. 

We show (A). The almost sure convergence holds because qei-) is a measurable function of lagged 
versions of r]{s,t) for s G Sm, t G Tt- Proposition 3 of Davis et al. [H] implies a strong law of large 
numbers. What remains to show is that the convergence is uniform on the compact parameter space 
0*. This can be done by carefully following the lines of the proof of Theorem 1 of Davis et al. [ 8 ], 
adapting it to the spatially anisotropic setting. For details we refer to Buhl [ 6 ], Theorem 4.4. We find 
that there is a positive finite constant Ki, independent of 0, M and T, such that 


E[|log5e(r/(s(^\i(^)),r/(s(2)T(2^))|] < ATi, fW), (s(2), 4 )) g (3.11) 


and that E[sup 0 g 0 * \qe{'^, 1; 'i',p)\ ] < 00 . Theorem 2.7 of Straumann m implies that the convergence 
is uniform. 


Next we show (B). Using Proposition 3 of Davis et al. [ 8 ] and (3.11) we have that, uniformly on 0*, 

P I 

E E E ‘“S{w (-!(». f)Ms + h.t + um 

/le-Hr «=o ' ^ {s,t)£gM,Tih,u) 


E 


as M, T —)• 00 . 


:[EE log {ge 1 ), f?(l + ^, 1 + u))} 

h^'Hr 

By Lemma and ( |3.11 ) it follows that, uniformly on 0*, 

p I 

E E le E log {»» (<l(o. *). 1(0 + h. * + «))} 


M'^r' 




{s,t)GgM,T(h,u) 


0 as M, T —)• 00 , 


Finally, we prove (C). Let 0 / 6*. For s G Sm and t G Tt, Jensen’s inequality yields 


E 


log 


ge {ri{s,t),g{s + h,t + u)) 
ge* ivis,t),g{s + h,t + u)) 


< log < E 


geigis,t),g{s + h,t + u)) 
ge* {g{s,t),g{s + hT + u))_ 


log I 


ge{vi,y 2 ) 


, W0*(2/i,y2)d(yi,y2)[ 

90 * (2/1,2/2) > 

log I j 90(2/1,2/2) d(yi,y2)} = 0 , 


{ 0 ,oo )2 


( 0 ,oo )2 


and it directly follows from (3.5) that PL(0) < PL(0*). As 9 T tlio identihability condition (3.10) 
yields (C). □ 
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ri 

r 2 

P 

identifiable parameters 

1 

0 

0 

Cl 

1 

1 

0 

Cl, C 2 

1 

1 

1 

Cl, C 2 , C 3 

> 1 

0 

0 

Cl, ai 

> 1 

> 1 

> 1 

Cl, ai, C 2 , Oi 2 , C 3 , as 


Table 1: Identifiable parameters for model (2.41 with d 


2 for some examples of maximum space-time lags (ri,r 2 ,p). 


Remark 1. There are combinations of maximum space-time lags that lead to non-identifiable param¬ 
eters, see Table However, Theorem still applies to all identifiable parameters (cf. Davis et al. [8], 
Remark 2). 

Next we prove asymptotic normality of the PMLE defined in (3.8). As in the proof of Theorem]^ 
we follow the lines of proof of Davis et al. [8], Section 5, adapting the arguments to the anisotropic 
setting. We start with some basic results needed throughout the remainder of the section. 

Lemma 2. Assume that all conditions of Theorem\^ are satisfied. Then for G and 

G [ 0 , 00 ), the following assertions hold componentwise: 

(1) The gradient of the bivariate log-density satisfies 

Ve logg0{r]{s^^\t^^'>), r]{s^‘^\t^‘^^)) 


E 


(2) The Hessian matrix of the bivariate log-density satisfies 


< 00 , 0 G 0 . 


E 


sup 

.flee* 




< 00 . 


Proof Assume identifiability of all parameters Cj, aj for j = 1,..., d -|- 1. For yi, 2/2 £ (0, 00 ) and for 
{h, u) G \ {0} lengthy but simple calculations of derivatives of yield 


V 0 logye(yi,y 2 ) = ^ VeS{h,u), 


d6{h, u) 


d6{h, u) 


dCi 


= \h 


^ log j = l,...d, 


and 


d6(h,u) I d6(h,u) ^ , , , 

V > ^ =\u\°‘d+l^ log|M|. 


dCd+i 


dad+i 


By compactness of the parameter space, as required in (3.9), we can bound those first partial derivatives 


as well as the second order partial derivatives from above and below. So it remains to show that for 
^(1) j s(2) g s and G T, 

aiog{y0(r/(s(^), r/(s0)^ t(2)))} |3n 


E^ 


and 


E, 


sup 

0ee* 


85{h, u) 

02 log{y0(r/(s(^), 72(5(2)^ t(2)))} 


< 00 


d‘^5{h,u) 


< 00 , 


where the function 6{h,u) can be treated as a constant since it is bounded away from 0 by (3.9). 
Hence, for the rest of the proof we refer to that of Davis et al. [8], Lemma 1. □ 






































For a central limit theorem we need certain mixing properties for a space-time setting (cf. Davis 
et al. [8], Section 5.1 and Huser and Davison m, Section 3.2). 

Definition 2 (Mixing coefficients and a-mixing). Let {r]{s,t) : s G Z'^,t G N} be a space-time process. 
Let d be some metric induced by a norm on For Ai, A 2 C x N let 

(i(Ai, A 2 ) := mf{d{{s^^\t^^'^), G Ai, {s^^\t^^'^) G A 2 }. 

(1) For k,i,n>0 the mixing coefficients are defined as 

ak/{n) :=snp{|P(Ai n A 2 ) -P(Ai)P(A 2 )| : 

Al G A 2 G F\ 2 , |Ai| < k, IA 2 I < i, (i(Ai, A 2 ) > n}, (3.12) 

where Fa. = cr{r]{s, t) : (s, t) G Aj) for i = 1,2. 

(2) {r]{s,t) : s G Z^,t G N} is called a-mixing if for all k,^>D, 

cik,e{n) —)• 0, n —>■ 00 . 


Recall from Eq. (2.2) that for {h,u) G with 5 as in (2.4) the tail dependence coefficient of 

the Brown-Resnick process is given by 


X{h,u) = 21 - 4 * 


\j 2 -I- • • • -I- -I- C'rf+i|t(|“'^+i] 


Corollary 2.2 of Dombry and Eyi-Minko |12) links the a-mixing coefficients with the tail dependence 
coefficients, and we will use this for the next result. 


Proposition 1. Let {r/(s,t) : s G G [0, 00 )} be the Brown-Resnick process (2.1) with dependence 
function 6 given by (|2.4[). Then the process {rj{s,t) : s G Z'^,t G N} is a-mixing, where the mixing 


coefficients in (3.12) satisfy for Hr as in (3.3) 


(1) n‘^ak,i{n) <oofork + l< ^{\Hr\ + l)(p+ 1), 

n=l 


(2) a(|^,|+l)(p+l),oo(?^) = o{n asn^oo, 

00 ^ 

(3) E '^(|'Hr| + l)(p+l),(|'Hr| + l)(p+l) ^ ^ 

n=l 

Proof. Note that for {h,u) G by the equivalence of norms, for some positive constant L, 


d{{h,u),{0,0)) < ^max{|/ii|,..., |/id|, |ri|} 

Therefore, for n G N, presuming d{{h, u), (0, 0)) > n results in max{|/ii|, ... ,\hd\, |rt|} > Ln, so that 
by Corollary 2.2 and Eq. (3) of Dombry and Eyi-Minko [12] we get 

C(k,e{'n) < 2k£ sup x{h,u) < 2ki sup x{h,u), (3.13) 

d{{h,u) ,{0,0))>n max{|/ii|,...,|/i(f|,|'u|}>Ln 

ak,oo{n) < 2k ^ x{h,u) < 2k ^ x{h,u). (3.14) 
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In the following we use the notation ||(h,,u)||oo := max{|/ii|,..., \hd\, |m|} for {h,u) G x N. Using 
1 — $(x) < exp{—for X > 0 and Eq. (2.2) and (3.13), we find for all k,i>0, 


Ok,i{n) < 4:k£ sup 1 — d>( 

\\{h,u)\\oo>Ln \ 


6{h, u) 


6{h, u) 


< 4k£ sup exp \ — 

\\{h,u)\\^>Ln 


= 4k£ sup exp<^--[Uil/iil"! + ...+ + C'rf+i|u|“‘*+i] 

\\(h,u)\\^>Ln I 4 

<Ak£ sup exp|-^min{|Ci|,...,|Cd+i|}||(h,u)||^“^"^’’''’"‘*+^^ 

\\(h,u)\\oo>Ln I 4 

< 4fc£exp |-J min{|Ci|,..., 




(3.15) 


—)■ 0 as n —)■ oo. 
This implies a-mixing. 


By similar arguments we obtain by (3.14) for all A; > 0, 

®fe,oo(^) — 

4A: ^ exp|-^min{|C'i|,...,|C'd+i|}||(/j.,'u) 

\\{h,u)\\ao>Ln 

We use the above bounds to prove assertions (l)-(3). 


OO 


(3.16) 


(1) For k + £ < 4{\1-Lr \ + l)(p + 1) we have by (3.15), 


OO OO ^ 

Y,n^au,t{n) < 4H J^n'^exp i -- min{|Ci|,..., |Crf+i|}(Ln)““^“i’-'“‘^+i> \ 

n=l n=l ^ ^ 


< OO. 


(2) First note that the number of grid points {h,u) G M^+i with ||(li,u)||oo = i for i G N equals 
(i + 1)'^^^ — and is therefore of order We use (3.16) and a more precise estimate than in 

part (1) to obtain for sufficiently large n 

^ *^(|'Ht'|+1)(p+1),oo(^) 

< 4n^+^(|Hr.| + l)(j>+ 1) 


\\{h,u)\\oo>Ln 


exp<^ --min{|C'i|,...,|C'rf+i|}||(h,M) 


I mm 
oo 




< 


OO .. 

Ksn<^^\\nr\ + l){p + l) Y *''exp{--min{Ci,...,Crf+i}i“-^“^’-’“'*+^>} 


i=[Ln\ 


0 as n —)■ oo. 


where is a positive constant. Convergence to 0 follows using the integral test for power series 
convergence and LemmaEq. (A.l). 

(3) We hnd, using again (3.15), 

OO 

Z]’^'^“(|w-i+i)(p+i),(iwd+i)(P+i)H" 

n=l 


10 












<(4[(|H.| + l)(p + l)]")^ 

OO - 

. ^ n'" exp { - - min{Ci,..., } 

n=l 


<OC 


as in ( 1 ). □ 

Because of Lemmaand Propositionthe following central limit theorem of Bolthausen [3] holds. 
Corollary 1. Consider the process {'V 0 qg*{s,t]r,p) : s G Z‘^,t G N}. Then 

1 ^ 

V 0 g' 0 *(s, t; r,p) ^ AA(0, Ei) as M, T —>■ oo, 

^'^y'^seSM t=i 

where 

OO OO OO 

E - E E Cov[V 0 g 0 *(l, l;r,p), V 0 g 0 *(si,.. .,Sd,t;r,p)] . (3.17) 

Si= —OO Sd = — 00 t=l 

Now we formulate the main result of this section. 

Theorem 2 (Asymptotic normality for large M and T). Assume the same conditions as in Theorem 
QJ Then 


y/M<^T{d - e*) ^ M{0, El) asM,T^ oo, 


where Ei := ^Ei(Fj^ ) with Ei given in (3.17) and 

Fi := E [-Vlq 0 *{l,l;r,p)] . 


(3.18) 


Proof. A Taylor expansion of the score function V around the true parameter vector 6 * 
yields for some 9 G [9,9*] : 

0 = VgPL^^’'^\9) = VePL^^''^\9*) + VIPL^^'^\9){9 - 9*). 


Therefore, 


M^Vf(9 -9*) = - 


M^T 


V: 


!p^(M,T)(0)y 


-1/ 1 


M 2 y/T 


VePL^^’^H9*)^ 




-1 


s^Sm ^—1 
T 


i=l 

= :-(/i-l2)-^(Ji-J 2 ). 

Note the following: 

• Corollary [^implies that Ji ^ A/’(0,Ei) as M,T —)■ 00 . 


E i2'^eq0<s,f,r,p) - ^—Ven^^’^\9*)) 
^ ^ MsVr ^ 
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Using representation (3.6) of the boundary term and Lemma we find 

p 


\\J2\\ = 


1 


< y^ 2 ' 


M2\/r 


EE E Ve \og{ge*{ri{s, t), r/(s + h,t + u))} 


mIVt 


1 


hS.S y/\QM,T{h,u)\ 

< y/K 2 {—= H- —) 

WM y/T 

P I 

hS.S V\SM,T{h,u)\ 


^ V0log{g0y{r]{s,t),r]{s + h,t + u))} 

{s,t)egM,T{h,u) 


^ V0log{g0y{g{s,t),g{s + h,t + u))} 

{s,t)egM,T(h,u) 


In the same way as done in Corollaryfor the process {Veg^* (s, t; r,p) : s G t G N}, we can 
apply Bolthausen’s central limit theorem to the processes {Ve logj^gj* {g{s, t),g{s + h,t + u))} : 
s G Z^, t G N} for G Hr, u € {0,... ,p}. We conclude that 


EE 


hGl-ir u=0 


^ \/\GM,Tih,u)\ 


^e'^og{g0*{g{s,t),vis + h,t + u))} 


(s,t)&gM,T{h,u) 


converges weakly to a normal distribution as M,T —)• oo, and it follows that J 2 —)• 0 as M,T —)■ 
00 . 

As {g{s,t) ■ s G Z‘^,t G N} is a-mixing, the process 

{'^0Qe{s,t-,r,p) : s G Z'^,t G N} 

is a-mixing as a set of measurable functions of mixing lagged processes. Furthermore, as 0 G 
[0,0*] and 0 is strongly consistent, we have that Ii —Fi as M,T —)> 00 . The convergence is 
uniform on 0* by Lemma which implies that 


E 


sup \Vlq0{l,l;r,p)\ 
0e0* 


< 00 . 


Concerning I 2 , the law of large numbers applied to 

I Vg \og{g 0 {g{s, t),g{s + h,t + ^i))} : s G Z"^, t G n| 

results in the fact that, in the same way as in part (B) of the proof of Theorem]^ I 2 0 as 
M, T —)> 00 . 


Finally, summarising these results, Slutzky’s Lemma yields (3.18). 


□ 
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3.2 Fixed spatial domain and increasing temporal domain 

As before we compute the PMLE based on observations on the area Sm x 7t, but now we consider 
M fixed, whereas T tends to infinity. 

We define the temporal a-mixing coefficients (cf. Ibragimov and Linnik |21j . Definition 17.2.1 or 
Bradley Hj, Definition 1.6). 

Definition 3 (Temporal mixing coefficients and temporal a-mixing). Let {r]{s,t) : s G SM,i £ hi} be 
a space-time process. Consider the metric d{-) of Definition^ 

(1) Let C N. For n > 0 the temporal a-mixing coefficients are defined as 

a{n) :=sup{|P(Ai n A 2 ) - P{Ai)P{A 2 )\ : 

Ai E , A 2 E P^^^^.j-( 2 ),d{SM X T^^\Sm X T^^^) > n}, (3.19) 

where PsmxP^) = t) € Sm x T^) for f = 1, 2. 

(2) {r]{s,t) : s E SmP £ N} is called temporally a-mixing, if 

a{n) —)• 0, n —>■ 00 . (3.20) 


Proposition 2. Let {rj{s,t) : s E E [0, 00 )} be the Brown-Resnick process (2.1) with dependence 
function 6 given by (2.4). Then the process {r]{s,t) : s E SmP E N} is temporally a-mixing, where the 
mixing coefficients (3.19) satisfy 


n=l 


n) 3 < oc. 


(3.21) 


Proof We use Eq. (3) and Corollary 2.2 of Dombry and Eyi-Minko |12] and (2.2) to obtain for n E N 


a{n) 

<2 sup EE 

d{SMxTW,SMxT^‘^'l)>n (^( 1 ) j(i)) ( 3 ( 2 ) j( 2 )) 

= 4 sup 

(i(5MXr(l),5AfXr(2))>n (s(2)_j(2)) 

(1 - $(^2 + • • • + Cfis^d^ - sf 1“^ + Cd+l\tm - t(2)|a.+i]) 

<4m 2'^ sup V 

d(5Mxr(i),5Mxr(2))>n (,(i) ,( 2 )) V ''V ^ ^ J 

er(i) xr(2) 

< 4M2'^ sup V exp I - 

<i(5Mxr(i),5Mxr(2))>n (,((^ 2 )) ^ ^ 

erW xr(2) 


where the last inequality follows from 1 — <l>(x) < exp{—for x > 0. We bound a(n) for large n 
further by 

a(n) < E E exp{-^Cd+i|t(i)-f(2)|«.+i|. 

t(i)e{ —oo,...,0} t(2)s{n,...,oo} 
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In the double sum a temporal lag u = > n appears exactly u — (n — 1) times. This yields 

OO - 

a{n) < 4M2‘^ ^(u - (n - 1)) exp | - ^Cd+iu '^‘^+^} 

u=n 

OO - 

< E u exp I - 


Convergence of the series (3.21) now follows by the integral test and Lemma 


□ 


In the following we show that strong consistency of the PMLE also holds, if the spatial domain 
remains fixed. 


Theorem 3 (Strong consistency for fixed M and large T). Assume the same conditions as in Theo¬ 
rem [7] restricted to the fixed space Sm ■ Then the PMLE 

0^ ’ ^ = argmax 


is strongly consistent, that is, 


'{M,T) 


as T ^ OO. 


Proof. For 0 G 0* and t GN, set 

p 

Q^{t]r,p) : = E E E l{(h,n)^{o,o)} log {ge (??(s, t),r]{s + h,t + u))} . 

s^Sm ■ii.=0 

s+hGS]^^ t+u<T 

Then 

PL^^^^\e) = Y,deif,r,p). 

t=i 

Following carefully the lines of the proof of Theorem the following conditions hold for fixed spatial 
domain: 

I T 

(A) — Yh Q.^{t',r,p) ^ PL^{6) := K[{q^[I] r,p)] as T —)• oo uniformly on the compact parameter 
T t=i 

space 0*. The main argument is that {■) is a function of temporally mixing lagged processes, 
then we apply again Theorem 2.7 of Straumann Ian. 

(B) The limit function PL^{6) is uniquely maximised at the true parameter vector 9* G 0*. 


□ 


Now we formulate the main result of this section. 

Theorem 4 (Asymptotic normality for fixed M and large T). Assume the same conditions as in 
Theorem restricted to the fixed space Sm ■ Then 

VT(e - 0*) ^ AA(0, Sa) asT^oo, (3.22) 

where Sa := with 

Fa:=E[-Vi(?^(l;r,p)] 
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and 


T.2 :=Yar[Veq^*il-,r,p)] + 2'^Cov[Vgq^*{l;r,p),V0q^*{t;r,p)]. 


t=2 

Proof. By its definition as a function of lagged temporally mixing processes, 
{'Vgq^{t;r,p))t^f^ is also temporally a-mixing with coefficients a'{n) = a{n — p). Furthermore, 

E [Ve log {ge* {r]{0, 0), ??(h, u))}] = 0, {h, u) £ 

because Lemma implies regularity conditions of the pairwise log-likelihood ( |3.4[ ) allowing to inter¬ 
change differentiation and integration. Now note that Lemma and Proposition imply that 

• KllVgq^{t; r,p)\^] < oo for f G N and every maximum spatial lag r and time lag p, and that 

OO ^ 

• Y1 Win)\3 < oo. 

n=l 

Therefore, the conditions of Theorem 18.5.3 of Ibragimov and Linnik m ( see also Bradley [3], Theo¬ 
rem 10.7) are satisfied and we conclude that 

T 


t=l 


r,p) ^ AA(0, S 2 ) as T —>■ 00 . 


(3.23) 


Taylor expansion of the score function VgPL^^''^\0) around the true parameter vector 6* yields for 
some 0 £[6,0*] ■. 


0 = VgPL^^’'^'>(e) = VgPL^^'^\e*) + vlPL^^'^\e)(e - e*). 


Therefore, 


\/r(0 -e*) = - (0)) ^ (0*)) 


t=i 


Vt 

^^J2'^oq0*{t-,r,p)) =: -T^J. 

t=i 


Note the following: 


(3.23) implies that J ^ AA(0,S2) as T —)• 00 . 

Uniform convergence holds because of Lemma which implies that componentwise 


E 


sup \Ylq^{l-,r,p)\ 
.6»e0* 


< 00 . 


By temporal a-mixing, since 6 £ [0,6*], and 6 is strongly consistent, we have I —F 2 as 
T —)• 00 . 


Finally, summarising those results, Slutzky’s Lemma yields (3.22). 


□ 


Throughout this section we have proved asymptotic properties of the parameter estimates of model 


(2.4) by classical results for ML estimators in combination with a spatio-temporal central limit theo¬ 
rem. Such results can also be applied to other models like geometrically anisotropic models, provided 
the required rates for a-mixing hold. 
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4 Test for spatial isotropy 


We use the results of Section to formulate statistical tests for spatial isotropy versus anisotropy 
based on the model (2.4), 


6{h,u) = Y,Cj\hjr^ +Cd+i\u\ 

i=i 


“d+i 


for spatial lags {h, u) = {hi,..., hd, u) G We derive the necessary results for d = 2. Generalisa¬ 

tions to higher dimensions are possible, but notationally much more involved. Again we consider the 
two cases of an increasing and fixed spatial domain. 


Dne to the structure of model (2.4) a test for isotropy versus anisotropy is a test of 
Hq : {Cl = C 2 and ai = 02 } versus Hi : {Ci 7 ^ C 2 or ai 7 ^ 02 }. 


(4.1) 


4.1 Increasing spatial domain 

From Theoremwe know that, under suitable regularity conditions, the PMLE 

6 = (Cl, C 2 , C 3 , «!, 02 , 03 ) 


is asymptotically normal; more precisely, for spatial observations on a regular grid and for T 
equidistant time points we have 


mVt 


(ci- cA 

C 2 -C 2 

C 3 -C 3 

Q-l — C(l 

a2 — 0.2 


ySs - a^j 


V 


N{Q, Si) as M, T —)• 00 , 


(4.2) 


where Si G is the asymptotic covariance matrix given in Theorem 

Our test is based on the spatial parameters only. Moreover, we test the two equalities in Hq 
separately and use Bonferroni’s inequality to solve the multiple test problem. 


Lemma 3. Assume the conditions of Theorem^ Setting Ai := (—1,1, 0, 0, 0,0) and A 2 := (0, 0, 0, —1,1, 0), 
we have that, as M, T —>• 00 , 


M Vf{{d2 - Cl) - (C 2 - Cl)) A AA( 0 , AitiAj), (4.3) 

M Vf{{a 2 - Si) - {02 - oi)) A J\f{0, A 2 tiAj). (4.4) 


Proof. We obtain the left hand side of (4.3) and (4.4) by multiplying Ai and A 2 to (4.2), respectively. 
This yields the limits on the right hand side by the continuous mapping theorem. □ 


We define 


6c := (C 2 — Cl), 6c := (C 2 — Ci), 6 a := {02 — ai), 6 a := (S 2 — Si). 
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Then the multiple test problem (4.1) becomes 


Ho,I : {9c = 0 } versus Hi,i : {Oq / 0 } 
Ho ,2 : {9a = 0 } versus Hi ,2 : { 6 »„ / 0 }. 


(4.5) 

(4.6) 


Since the variances in (4.3) and (4.4) are not known explicitly, we find the rejection areas of the two 

I, Chapter 5. Their main Assumption 5.3.1, the 


tests by subsampling as suggested in Politis et al. 
existence of a weak limit law of the estimates, is satisfied by Lemma 

We formulate the subsampling procedure in the terminology of the space-time process { t /( s , t) : s € 
SmH G 7t}- We choose space-time block lengths 

h = {hi,b2-,bo) > (1,1,1) and the degree of overlap e = (61,62,63) < {M, M,T). The blocks are 
indexed by i = {ii,i 2 ,i 3 ) G with ij < qj for qj := [J -|- 1, j = 1,2 and qo := \ ^--^\ + 1- This 
results in a total number oi q = qiq 2 Q 3 blocks, which we summarise in the set 


Ei,b,e = {(si,S 2 ,t) GSm xTt : {ij - l)ej -b 1 < sj < {ij - l)ej +bj,j = 1 , 2 , 
{is - 1)63 -b 1 < f < (^3 - 1)63 + 63}- 


Now we estimate 9c and 9a based on all observations in a block, hence getting q different estimates, 
which we denote by 9c,b,i and 9a,b,i- 

In order to find rejection areas for the isotropy test, we will use Lemma and take care of the 
unknown variance in the normal limit by a subsampling result. 

Theorem 5. Denote by tm,t ■= My/T and the square roots of the number of observa¬ 

tions in total and in eaeh block, respectively. Assume that the conditions of Theorem hold and, as 
M, T —>■ 00 , 

(i) —)• 00 for i = 1, 2, 3, such that hi = o{M) for i = 1,2, and 63 = o{T) (hence, tij/tm,t —t OJ, 

(a) e does not depend on M or T. 

In the following 9 stands for either 9c or 9a. Define the empirical distribution function 

^ Ql 92 93 

^ * 1=1 *2 = 1 *3 = 1 

and the empirical quantile funetion 

'^ 6 , 0(1 “ ^) ■= I® ^ ® 1 ’ /3e(o,i). (4.8) 

Then the following statements hold for M, T —)■ oo.' 

(1) Denote by 4>ct(-) the distribution function of a mean 0 normal random variable Z with variance 

2 AiIAAj, in ease of 9c, 

I A 2 T 11 A 2 , in ease of 9a, 

and recall that 24>(j(-) — 1 is the distribution function of |Z|. Then 

Lbg{x) 24>ct(x) — 1, a: G M. 
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(2) Set Jn{x) := IP(tm,t|^ — < a:) for x G M, then 


sup 




(3) For (3 G (0,1), 


TM,t\ 0 - d\ < Cf, 0 (l - 


0 . 


1-/3. 


(4.9) 


Proof. We apply Corollary 5.3.1 of Politis et al. [29]. Their main Assumption 5.3.1; i.e., the existence 
of a continuous limit distribution, is satisfied by Lemma Assumptions (i)-(ii) are also presumed by 
Politis et al. [29]. The required condition on the a-mixing coefficients is satisfied similarly as in the 
proof of Proposition by Lemma [^ and the result holds. □ 


From (4.9), we find rejection areas for the test statistics at confidence level /3 G (0,1) as 

(recall that 9 stands for either 9c or 9a) 

:= (- 00 , - /3)) U - /3), oo) = - P), c^ g{l - /3)]F 

Bonferroni’s inequality 

^(reject iLo,i or iLo, 2 ) < IP(reject iLo,i) + IP(reject iLo, 2 ) < 2/3, 
applies and solves the multiple test problem. 


4.2 Fixed spatial domain 

First note that an analogue of Lemma holds with rate \/T instead of My/T and with the asymptotic 
covariance matrix S 2 as given in Theorem 

The subsampling statement corresponding to Theorem then reads as follows. 

Theorem 6. Denote by tt ■= VT and = y/b^ the square roots of the number of time points of 
observations in total and in each block, respectively. Assume that the conditions of Theorem are 
satisfied and that Lemma^holds for T —)• 00 with rate y/T instead of My/T and with the asymptotic 
covariance matrix S 2 as given in Theorem\^ Assume further that as T ^ 00 , 

(i) 63 —)• 00 such that 63 = o{T) (hence, > Oj, 

(ii) e does not depend on T, 


(Hi) bi, 62 —)• M. 

Let b = ( 61 , 62 ,^ 3 ); Xb = y/Fjxibz and tm,t = My/T. With Si as in Theorem^ replaced by M^S 2 , 
conclusions (a), (b), and (c) of Theorem^ remain true as T tends to infinity. 

Proof. We apply Corollary 5.3.2 of Politis et al. [221. The required temporal mixing condition is 
satisfied similarly as in the proof of Proposition by Lemma [^ □ 


Remark 2. We can in practice apply the same procedure of subsampling as in Section 4.1 This is 
justified by the fact that Tb^jxT —^ 0 implies that Tb/TM,T —>■ 0 as T —)• 00 under conditions (|^-(iii) 


of Theorem]^ In particular, the rejection area for tt9 (where again 9 stands for either 9c or 9a) is 
found as 
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Figure 1: Rainfall observation area in Florida 


5 Data analysis 


We fit the Brown-Resnick space-time process (2.1) with dependence structure given by the model (2.4) 


to radar rainfall data, which were provided by the Southwest Florida Water Management District 
(SWFWMD). The data used for the analysis are rainfall measurements on a square of 120kmx 120km 
in Florida (see Figure over the years 1999-2004. The raw data consist of measurements in inches on 
a regular grid in space every two kilometres and every 15 minutes. Since there exist wet seasons and 
dry seasons with almost no rain we consider only the wet season June-September. Moreover, the area 
is basically flat with predominant easterly winds due to its closeness to the equator and, therefore, 


existing trade winds. Hence, (2.4) with parameters that possibly differ along both spatial axes fits well 


without introducing a rotation matrix. 


5.1 Data transformation and marginal modelling 

We carry out a block-maxima method in space and time as follows: We calculate cumulated hourly 
rainfall by adding up four consecutive measurements. Then we take block-maxima over 24 consecutive 
hours and over lOkmxlOkm areas; i.e., the daily maxima over 25 locations, resulting in a 12 x 12 
grid in space for all 6 x 122 days of the wet seasons giving a time series of dimension 12 x 12 and 
of length 732. Taking smaller areas than lOkmxlOkm squares or a higher temporal resolution (e.g. 
12 -hour-maxima) results in observations that are not max-stable and the max-stability test described 
in Section would reject. 

By removing possible seasonal effects, we transform the data to stationarity. We obtain the obser¬ 
vations 


{f)((si, S2),t) : si, S2 = 1,..., 12, f = 1,..., 732} . 


(5.1) 
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Taking daily maxima removes for every location most of the dependence in the time series. This 
implies that marginal parameter estimates found by maximum likelihood estimation are consistent 
and asymptotically normal. 

To give some details: for each fixed location (si,S 2 ), we fit a univariate generalised extreme value 
distribution (cf. Embrechts et al. [H], Definition 3.4.1) to the associated time series. The estimated 
shape parameters are all sufficiently close to 0 to motivate a Gumbel distribution as appropriate model. 
We therefore fit a Gumbel distribution A^^o-{x) = exp{—e with parameters fi = ^(si,S 2 ) £ M 

and a = cr(si, S2) > 0 and obtain estimates = /u(si, S 2 ) and a = a{si, S2). 


Depending on different statistical questions and methods, we transform (5.1) either to standard 
Gumbel or standard Frechet margins. In the first case we set 


m((si,'S 2 )T) : = 


fji{si,S2),t) - /i 


t = 1, 


a 


,732, 


and in the latter case, with j denoting the Gumbel distribution with estimated parameters, 

1 


V2iisi,S2),t) := - 


log{A^,5(f?((si,S2)T))}’ 


t = 1,...,732. 


(5.2) 


(5.3) 


We assess the goodness of the marginal fits by qq-plots of the observations (5.2) versus the standard 
Gumbel quantiles for every spatial location. Figure depicts the qq-plots at four exemplary spatial 
locations (1,1), (6,8), (9,4) and (11,10). Gonfidence bounds are based on the Kolmogorov-Smirnov 
statistic (cf. Doksum and Sievers [llj . Theorem 1 and Remark 1). All graphs show a reasonably good 
fit. 


In the following data analysis we regard (5.3) as realisations of the space-time Brown-Resnick 
process (2.1) with dependence structure 5 as in (2.4): 


d{huh2,u) = ci\hr +C2\h2r+C3\ur, ( 5 . 4 ) 

with hi = /i 2 = ^ 2 ^^ — u = — t^'^'l, for two spatial locations = (s^^\s 2 ^^) and 

and two time points t^^^ and t^‘^'1. 


5.2 Testing for max-stability in the data 

We first want to check if the block-maxima data originate from a max-stable process. A diagnostic 
tool is based on a multivariate Gumbel model (cf. Gabda et al. [E]), and we explain first the method 
in general. We assume a space-time model of a general spatial dimension d G N. As before, we denote 
the regular grid of space-time observations by 

5MxrT = {i,...,M}''x{i,...,r}. 


We define a hypothesis test based on the standard Gumbel transformed space-time observations (5.2) 
by 


Ho : {r]i{s,t) : {s,t) G x [0,cx))} is max-stable. 
^We use the R-package extRemes (Gilleland and Katz |17|i. 


(5.5) 
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Location (1,1) 


Location (6,8) 




Gumbel 


Gumbel 


Location (9,4) 



Gumbel 


Location (11,10) 



Gumbel 


Figure 2: qq-plots of the Gumbel transformed time series values versus the standard Gumbel distribution for four 
locations: (1,1) (top left), (6,8) (top right), (9,4) (bottom left) and (11,10) (bottom right). Dashed blue lines mark 95% 
confidence bounds. Solid red lines correspond to no deviation. 


Under Hq all finite-dimensional margins are max-stable; particularly, for every D C Sm x Tt, the 
multivariate distribution function of {r/i(s,t) : (s, t) G D} is given by 


Gniyi, ■ ■ ■,y\D\) = exp{-VD(e^^ . • • ,e^i^l)}, (t/i,.. .,y\D\) G 


where Vd is the exponent measure from (2.3). Since Vd is homogeneous of order -1, the random 
variable 


rji) := max{? 7 i(s,t) : {s,t) G D} 
has univariate Gumbel distribution function 


< y) = GdC?/, ... ,y) = exp{-e ^yD(l,... ,1)} = e ® yeR; (5.6) 


i.e., fj-D '■= log Vd( 1, •.., 1) is the location parameter and, since 1 < Vd(1,...,1) < \D\, we have 
0 < fJ-D ^ log|U*|. These considerations can be used to construct a graphical test for max-stability: 
First, choose different subsets D with the same fixed cardinality. Then extract several independent 
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Figure 3: R independent realisations oi rjo for different subsets D of the space-time observation area. 


realisations of the random variables !]£> from the data and test by means of a qq-plot, if they follow a 
Gumbel distribution. 


We apply this test to the standardized Gumbel transformed data (5.2). As indicated above, taking 


daily maxima removes for every location most of the dependence in the time series. For this test we 
want to take every precaution to make sure that we work indeed with independent data. Preliminary 
tests show that spatial observations, which are a small number of B 2 days apart (to be specified below), 
show only very little time-dependence. 

Consequently, we define time blocks of size Bi of spatial observations, which are in turn separated 
by time blocks of size B 2 as 


Sm X r® = {1 ,..., M}2 X {(i - l)(Si + B 2 ) + t : t = 1 ,..., Bi}, 


(5.7) 


for i = 1,..., i? = ■ The numbers Bi and B 2 need to be chosen in such a way that the blocks 

can be considered as independent. This results in R independent time blocks of length Bi of spatial 
data and thus in R independent realisations of !]£> for every D C Sm x {1,..., i?i}. The procedure is 
illustrated in Figure]^ 

We use these i.i.d. realisations to estimate for every D by maximum likelihood estimation 
restricted to [0,log|T)|]. Since the MLE of the location parameter of a Gumbel distribution is not 
unbiased (cf. Johnson et al. |22j . Section 9.6), we perform a bias correction. 

For the diagnostic we take X G N and consider subsets D with cardinality \D\ = K. As the 
total number ( ^^ ) of those subsets is in most cases intractably large, we randomly choose m := 
min{i?, )} subsets and obtain in total N = m-R subsets, which we denote by for j = 1,..., m 
and i = 1,..., R. For every j = 1,..., m we estimate by MLE based on the i.i.d. random variables 
■= Vr)(i ), * = 1, • • •, J?- Then we perform qq-plots of 

(1) (m) 

RDrn 
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qq-plot for K=2 


qq-plot for K=3 




Gumbel 


Gumbel 


qq-plot for K=4 



Gumbel 


qq-plot for K=5 



Gumbel 


Figure 4: qq-plots of theoretical standard Gumbel quantiles versus the empirical quantiles (black dots). The latter 
correspond to the empirical distribution of maxima taken over groups of cardinality K. Dashed blue lines mark 95% 
pointwise confidence bounds obtained by block bootstrap. Solid red lines correspond to no deviation. 

versus the standard Gumbel distribution. As a measure of variability of the estimates, non-parametric 
block bootstrap methods (cf. Politis and Romano [28j, Section 3.2) are applied to obtain 95% pointwise 
confidence bounds. Using bootstrap methods, we preserve the dependence between different subsets 
D in the confidence intervals. Under Hq, the bisecting line should lie within these confidence bounds. 

The Florida daily rainfall maxima show only little temporal dependence beyond one day. Hence 
we choose Bi = 2 and B 2 = 1, which yields R = = 244 mutually independent time blocks of 

spatial data. We perform the described procedure for K = 2,3, 4, 5, which entails m = R = 244. Thus 
we obtain a total number of = 244^ = 59 536 subsets. The power of this diagnostic test increases 
with K (cf. Gabda et al. [E]) as it gets less likely to include sets of space-time points that are iF-wise 
independent. Figure shows the results for the different choices of K. The solid red bisecting lines 
lie inside the confidence bounds. Hence, there is no statistically significant evidence of the space-time 
process generating the data not to be max-stable. 
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5.3 Pairwise maximum likelihood estimation 


We apply the pairwise maximum likelihood estimation to the standard Frechet transformed data 
(5.3). The parameters to estimate are those of the function 5 in (5.4); i.e., C'i,C 2 ,C '3 G (0,oo) and 
«!, 02 , as G (0, 2]. 


In the definition of the pairwise log-likelihood function (3.4), the maximum spatial and temporal 


lags are specified by the numbers ri, r 2 and p, respectively. Immediately by model (5.4) for 6 , the 
parameters of the three different dimensions (space and time) are separated in the extremal setting. 
This has also been noticed in Davis et al. [ 8 ], where a simulation study in Section 7 for the isotropic 
model shows that estimating the spatial and temporal parameter pairs individually leads to very 
good results in terms of root-mean-square error and mean absolute error. Hence, for example for 
parameter estimates for Ci and ai, we can set the maximum lags corresponding to the remaining 
parameters equal to 0 (i.e., we set r 2 = p = 0). This means that we basically fit univariate models 
to the respective spatial and temporal parts of the dependence function (5.4). Hence, this separation 
simplifies the statistical estimation. However, proving asymptotic properties of the pairwise likelihood 
estimator in the special case of a univariate model would for instance still involve showing the required 
mixing conditions and thus not remove much of the complexity. 

Furthermore, we know that we should not include too many lags in space or time into the likelihood, 
since independence effects can introduce a bias in the estimates, see for example Nott and Ryden [25] . 
Section 2.1, or Huser and Davison m, Section 4. On the other hand, an empirical analysis showed 
that extremal spatial dependence of the Florida daily rainfall maxima ranges up to lag 4 and extremal 
temporal dependence does not last more than one or two days, cf. Figure 7.2.6 in Steinkohl m- Hence, 
we perform the PMLE for maximum spatial and temporal lags up to 4 and 2, respectively, thus also 
assuring identifiability of all parameters according to Table The results are summarised in Table 
Setting ri, r 2 or p equal to 1 results in non-identifiability of the corresponding parameters ai, 02 or 
as, respectively; cf. TableTherefore, they are not shown in Table 

The combination of a rather large estimate for C 3 and a rather small estimate for as indicates 
that there is only little extremal temporal dependence, see Steinkohl [30], Section 7.2. Asymptotic 
95%-confidence intervals are based on asymptotic normality of the parameter estimates and estimated 
using subsampling methods (cf. Section . 


5.4 Isotropic versus anisotropic model 


Using the results of Section]^ we want to apply the test (4.1) for spatial isotropy to the hypothesis 


Hq : {Cl = C 2 and ai = 02 } versus Hi : {Ci ^ C 2 or 01 / 02 }. 


For the block maxima of the precipitation data we have d = 2, M = 12 and T = 732. This 
corresponds to the situation of a fixed spatial domain with tt = \/732. 

We use the spatial PMLEs based on maximum lags 2-4, which can be read off from Table We 
obtain the rejection areas from Theorem]^ We choose hi = 62 = 5, thus ensuring that the full range 
of spatial dependence is contained in the subsamples and simultaneously achieving that their number 
is large. Concerning the number of time points in each subsample, we take 63 = 600. Here we choose a 
large number to ensure that Theorem]^ where T —)• 00 , is applicable. This results in = -v /63 = \/600. 
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max. lags 

a 

Si 

( 2 , 0 , 0 ) 

0.6287 

[0.5928, 0.6646] 

0.9437 

[0.9065,0.9808] 

(3,0,0) 

0.6358 

[0.5989, 0.6728] 

0.8599 

[0.8189,0.9009] 

(4,0,0) 

0.6438 

[0.6051,0.6825] 

0.8107 

[0.7690,0.8525] 

( 0 , 2 , 0 ) 

0.7271 

[0.6492, 0.8050] 

0.9517 

[0.8715,1.0320] 

(0,3,0) 

0.7370 

[0.6586,0.8154] 

0.8521 

[0.7737,0.9305] 

(0,4,0) 

0.7476 

[0.6677, 0.8275] 

0.7931 

[0.7039,0.8822] 

( 0 , 0 , 2 ) 

4.8378 

[4.4282, 5.2474] 

0.1981 

[0.0177,0.3784] 


Table 2: Estimates of the parameter pairs {Ci,ai), ( 02 , 0 . 2 ) and (( 73 , 03 ) for different maximum spatial and temporal 
lags. Intervals below the point estimates are asymptotic 95%-confidence bounds based on subsampling. 


In order to obtain a large number of subsamples, we further choose ei = 62 = 63 = 1 as the degree of 
overlap. 


max. 

lag 

Tt 

C 2 -C 1 

tt{C2 - Cl) 

Rej^^'^ 

■'dc 

97.5%-CI 

for C 2 - Cl 

Reject 

Cl=C2 

2 

27.055 

0.098 

2.651 

[-2.400,2.400]= 

[0.010,0.187] 

yes 

3 

27.055 

0.101 

2.738 

[-2.392,2.392]= 

[0.013,0.190] 

yes 

4 

27.055 

0.104 

2.808 

[-2.393,2.393]= 

[0.015,0.192] 

yes 


Table 3: Test results for parameters (7i and (72. All values are rounded to three positions after decimal point. 


max. 

lag 

Tt 

02 ~ Ctl 

TT{oi2 — Ctl) 

RejC^'^ 

97.5%-CI 

for ( 0:2 — oi) 

Reject 

Ctl = OL 2 

2 

27.055 

0.008 

0.216 

[-2.162,2.162]= 

[-0.072,0.088] 

no 

3 

27.055 

-0.008 

-0.216 

[-2.130,2.130]= 

[-0.087,0.071] 

no 

4 

27.055 

-0.018 

-0.477 

[-2.342,2.342]= 

[-0.104,0.069] 

no 


Table 4: Test results for parameters oi and 02 . All values are rounded to three positions after decimal point. 


Tables and 1^ present the results of the two tests at individual confidence levels /3 = 2.5% giving 
a test for ( |4.1[ ) at a confidence level 2(3 = 5% by Bonferroni’s inequality. The differences (C 2 — Ci) 
and (02 — Si) can be obtained from Table 

Since we can reject the individual hypothesis that Ci = C 2 at a confidence level of 2.5%, we 
can reject the overall hypothesis Hq of (4.1) at a confidence level of 5% and conclude that our data 
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originate from a spatially anisotropic max-stable Brown-Resnick process. Further note the interesting 
fact that, although the asymptotic confidence interval for the difference C 2 — Ci does not include 0, 
the individual intervals for Ci and C 2 overlap, see Table This is due to the fact that the individual 
confidence bounds are estimated independently of each other, whereas the estimated bounds for the 
difference reflect how far the parameter estimates lie apart in one fixed particular (sub)sample. 


5.5 Model check 


Finally, having fitted the Brown-Resnick space-time model (2.1) to the precipitation data, we want to 


assess the quality of the fit. We take inspiration from Section 5.2 of Davison et al. [9] and compare 
maxima taken over subsets of the space-time precipitation data with simulated counterparts. 


Similarly as in Section 5.2, we consider subsets of the observations on a regular grid for L spatial 
locations and for time points 1,..., Ri, 






Ri) L}. 


We follow the procedure as in (5.7) to extract R independent realisations of {r/i(s, t) : (s, t) G D] from 


the standard Gumbel transformed space-time observations (5.2). This yields in turn R independent 


realisations of r]D = max{r 7 i(s,t) : (s,t) G R}, which we summarise in the ordered vector r/jata • = 

,..., Now we simulate a corresponding vector, denoted by rysim := (4?^ • ■ • > 4^4 To this 

end we need reliable Monte Carlo values as elements of i^sira- We obtain them by simulating empirical 

order statistics as follows. We simulate m ■ R independent copies of the Brown-Resnick space-time 

process on D with dependence structure 5 as in ( |2.4[ ) with the PMLEs from Tablewhere we take the 

estimates based on maximum lag 4 (for the spatial parameters) and 2 (for the temporal parameters), 

which are the maximum lags, where dependence is still present. We transform the univariate margins 

to standard Gumbel. This results in corresponding m-R independent simulations of and we consider 

^(i) 

them as m blocks of size R. We order the R values in each block and define as the mean of all 
simulated ith order statistics for i = 1,..., R, which gives rjsim ■= ■ ■ ■ 1 4^4 

The vectors r/data and rjsim are compared by qq-plots. If the fit is good, the points in the plots lie 
approximately on the bisecting line. Pointwise 95%-confidence bands are determined by the 2.5% and 


the 97.5% quantiles of the simulated order statistics. As in Section 5.3, we choose Ri = 2. The number 


of simulations is iV = m ■ R = 100 • 244 = 24400. Figure presents the results for four exemplary 
groups of locations. The plots reveal a good model fit. 

We carried out the simulations using the exact method recently suggested in Dombry et al. m, 
Sections 3.3 and 5.2. For an overview and comparison of different simulation methods for Brown- 
Resnick processes we refer to Leber 


5.6 Application: conditional probability fields 

Based on the fitted model, we want to answer questions like: Given there is extreme rain at some 
space-time reference point (s^, S 2 , t*) G {1,..., 12}^ x {1,..., 732}, what is the estimated probability 
of extreme rain at some prediction space-time point ( 5 ^, 53 , t^)? In other words, we want to estimate 
the probabilities 

> Z I v{{s*,S 2 ),t*) > z*) , (5.8) 
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GOF qq-plot for L=2 


GOF qq-plot for L=3 




GOF qq-plot for L=4 GOF qq-plot for L=5 




Figure 5: Goodness of fit qq-plots for different spatial locations and different L. Top left: L — 2-. (1,1) a^nd (1,2). Top 
right: L = 3: (1,1), (1,2) and (3,1). Bottom left: L = 4: (1,1), (1,2), (3,1) and (3,2). Bottom right: L = 5: (1,1), (1,2), 
(3,1), (3,2) and (2,1). PMLEs underlying the simulations are based on maximum spatial and temporal lags 4 and 2, 
respectively. Dashed blue lines mark 95% pointwise conhdence bounds. Solid red lines correspond to no deviation. 


where S 2 ),t) : si, S 2 = 1, • • ■, 12, t = 1,..., 732} are the stationary observations (5.1) and z and 

z* are prediction and reference rainfall levels, respectively. Denote by the Gumbel distribution 


with location and scale parameters fj. and a (cf. Section 5.1) and set JP := ^{ 3 ^, 82 ), := u(s^,S 2 ), 


Jl* := /2 (s},S2 ) and u* := a(s},S 2 ), which are the marginal Gumbel parameter estimates. Simple 


computations show that (5.8) can be estimated by 


1 - (z*) 

+ exp I - Gd 


1 


log{A^p^S:p(2:)} ’ log{A^*^5*(2;*)} 


)}). 


where Vo is the estimate of the exponent measure ( |3.2[ ) obtained by plugging in the PMLEs of the 
parameters of the dependence function 6 . Figure shows four predicted conditional probability fields 
for the reference points (1,1,1), (5, 6,1), (8,10,1) and (10, 7,1) and for high empirical rainfall levels 
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Figure 6: Predicted conditional probability fields based on daily maxima for reference space-time points (1,1,l)i (5,6,1), 
(8,10,1) and (10,7,1) and rainfall levels z = z* = 2.5 (clockwise from the top left to the bottom right). 

2 : = z* = 2.5. Because of the little temporal dependence in the daily maxima, we only consider equal 
time points for spatial predictions. 

Acknowledgements 

We take pleasure in thanking Anthony Davison and his group for an extremely pleasant and interesting 
time at EPFL Lausanne, and SB also at a Summer School in Leukerbad. SB also thanks Christina 
Steinkohl for her constant support, when writing his Master’s Thesis. Discussions with Richard Davis 
and Jenny Wadsworth are gratefully acknowledged. We also thank Chin Man Mok for providing the 
Florida rainfall data and acknowledge that the data are provided by the Southwest Florida Water 
Management District (SWFWMD). We thank the referees for fruitful comments and remarks. SB 
additionally acknowledges that he was supported by Deutsche Forschungsgemeinschaft (DFG) through 
the TUM International Graduate School of Science and Engineering (IGSSE). 

References 

[1] J. Beirlant, Y. Goegebeur, J. Segers, and J. Teugels. Statistics of Extremes, Theory and 
Applications. Wiley, Chichester, 2004. 


28 





























[2] J. Blanchet and A. Davison. Spatial modeling of extreme snow depth. Ann. Appl. Stat., 5(3): 


1699-1724, 2011. 

[3] E. Bolthausen. On the central limit theorem for stationary mixing random fields. Ann. Probab., 
10(4):1047~1050,1982. 

[4] R. Bradley. Introduction to Strong Mixing Conditions, volume I. Kendrick Press, Heber City, 
Utah, 2007. 

[5] B. Brown and S. Resnick. Extreme values of independent stochastic processes. J. Appl. Probab. , 
14(4);732-739, 1977. 

[6] S. Buhl. Modelling and Estimation of Extremes in Space and Time. Master’s thesis, Technis- 
che Universitat Miinchen, 2013. Available under https://niediatum.ub.tum.de/doc/1145694/ 
1145694.pdf, 

[7] R. Davis, C. Kliippelberg, and C. Steinkohl. Max-stable processes for modeling extremes observed 
in space and time. J. Korean Stat. Soc. , 42(3):399-414, 2013. 

[8] R. Davis, C. Kliippelberg, and C. Steinkohl. Statistical inference for max-stable processes in space 
and time. J. Roy. Stat. Soc. B , 75(5):791-819, 2013. 

[9] A. Davison, S. Padoan, and M. Ribatet. Statistical modelling of spatial extremes. Stat. Sci. , 27 
(2):161-186, 2012. 

[10] L. de Haan. A spectral representation for max-stable processes. Ann. Probab., 12(4);1194-1204, 
1984. 

[11] K. Doksum and G. Sievers. Plotting with confidence: graphical comparisons of two populations. 
Biometrika , 63(3):421-434, 1976. 

[12] C. Dombry and F. Eyi-Minko. Strong mixing properties of max-infinitely divisible random fields. 
Stoch. Proc. Appl. , 122(11):3790~3811, 2012. 

[13] C. Dombry, S. Engelke, and M. Oesting. Exact simulation of max-stable processes. Biometrika , 
103, 2016. To appear. 

[14] P. Embrechts, C. Kliippelberg, and T. Mikosch. Modelling Extremal Events. Springer, Berlin, 
1997. 

[15] S. Engelke, A. Malinowski, Z. Kabluchko, and M. Schlather. Estimation of Hiisler-Reiss distri¬ 
butions and Brown-Resnick processes. J. Roy. Stat. Soc. B , 77(l):239-265, 2015. 

[16] D. Gabda, R. Towe, J. Wadsworth, and J. Tawn. Discussion of ” Statistical Modelling of Spatial 
Extremes” by A.C. Davison, S.A. Padoan and M.Ribatet. Stat. Sci., 27(2):189-192, 2012. 

[17] E. Gilleland and R. W. Katz. New software to analyze how extremes change over time. Eos, 92 
(2):13-14, 2011. 


29 


















[18] E. Gine, M. G. Hahn, and P. Vatan. Max-infinitely divisible and max-stable sample continuous 
processes. Probab. Theory Rel. , 87:139-165, 1990. 

[19] R. Huser and A. Davison. Space-time modelling of extreme events. J. Roy. Stat. Soc. B , 76(2): 
439-461, 2014. 

[20] J. Hiisler and R.-D. Reiss. Maxima of normal random vectors: between independence and complete 
dependence. Stat. Probabil. Lett. , 7(4):283-286, 1989. 

[21] I. Ibragimov and Y. Linnik. Independent and Stationary Sequences of Random Variables. Gronin¬ 
gen: Wolters-Noordhoff, 1971. 

[22] N. Johnson, S. Kotz, and N. Balakrishnan. Gontinuous Univariate Distributions - Volume 2 . 
Wiley, New York, 2nd edition, 1995. 

[23] Z. Kabluchko, M. Schlather, and L. de Haan. Stationary max-stable fields associated to negative 
definite functions. Ann. Probab. , 37(5):2042-2065, 2009. 

[24] D. Leber. Gomparison of Simulation Methods of Brown-Resnick processes. Master’s thesis, 
Technische Universitat Miinchen, 2015. Available under https://mediatuni.ub.tum.de/doc/ 
1286600/1286600.pdf, 

[25] D. J. Nott and T. Ryden. Pairwise likelihood methods for inference in image models. Biometrika, 
86(3):661-676, 1999. 

[26] M. Oesting. Simulationsverfahren fiir Brown-Resnick-Prozesse. Diploma thesis, Georg-August- 
Universitiit zu Gottingen, 2009. Available under http://arxiv.org/abs/0911.4389v2, 

[27] S. Padoan, M. Ribatet, and S. Sisson. Likelihood-based inference for max-stable processes. E 
Am. Stat. Assoc. , 105(489):263-277, 2010. 

[28] D. N. Politis and J. P. Romano. Nonparametric resampling for homogeneous strong mixing 
random fields. J. Multivariate Anal., 47:301-328, 1993. 

[29] D. N. Politis, J. P. Romano, and M. Wolf. Subsampling. Springer, New York, 1999. 

[30] C. Steinkohl. Statistical Modelling of Extremes in Space and Time using Max-Stable Processes. 
Dissertation, Technische Universitat Miinchen, 2013. Available under https://mediatum.ub. 
turn.de/doc/1120541/1120541.pdf, 

[31] D. Straumann. Estimation in Conditionally Heteroscedastic Time Series Models. Lecture Notes 
in Statistics, No. 182, Springer, Berlin, 2004. 

[32] J. Wadsworth and J. Tawn. Efficient inference for spatial extreme value processes associated to 
log-gaussian random functions. Biometrika, 101(1):1-15, 2014. 

[33] A. Wald. Note on the consistency of the maximum likelihood estimate. Ann. Math. Stat., 20(4): 
595-601, 1949. 


30 



















A An auxiliary lemma 

Lemma 4. The following two bounds hold true for r > 1, a € (0,2] and C > 0; 


f u^'e du ~ 




2/ OO, 


/ L\J \ 

f ( f du j ^ dy < oo. 

1 ^ 


(A.l) 

(A.2) 


Proof. First note that integrals of the form drt are finite for every r > —1, a G (0, 2], and 

C > 0, since they are transformations of the gamma function r(x) = dt, which exists for 

positive X. We prove (A.l) by an application of rHopital’s rule: 


°° .,,r^-Cu° 


f °° u^e 
lim ^ 


du 


= lim 


-y^Q-Cy^ 


= lim 


y^oo (^-yr _g ^-^+^ y^-«) Q-Cy^- y^oo yT (l _ 


= 1 . 


In order to prove (A.2) first note that it follows from ( |A.1[ ) that for every e > 0 there exists yo = yo(e) 
such that for all y > yo, 


1 1 i 

r —Cu^ r-a+1 _C o 

u e duj <(l + e)(^—j y s e 3^ 


y 


Now we split the double integral of (A.2) up into 

yo oo ^ oo 

yT ^-CU°‘ ^ fij/ + 


u''e du) ^ dy =: Ii + 12 - 


For Ji we obtain 


1 y 


yo yo 


yo y 


yo oo 


A < y ( y du) ® dy + J (^J du) ® dy =: + lf\ 

I y I yo 

is obviously finite, and to bound we use (A.3), which yields 

1 \ 5 r-a+l 


Concerning I 2 , note that 


< (2/o-l)(l + e)(^)''yo " e s^o <c 00 . 


CXJ 

/2 < (1 + ^) (^) " y 2/"^^e"f dy. 


yo 


which is finite by finiteness of the gamma function. 


(A.3) 


□ 
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